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Abstract 

We propose a compressive sensing algorithm that exploits geometric properties of images to recover images 
of high quality from few measurements. The image reconstruction is done by iterating the two following steps: 1) 
estimation of normal vectors of the image level curves and 2) reconstruction of an image fitting the normal vectors, 
the compressed sensing measurements and the sparsity constraint. The proposed technique can naturally extend to 
non local operators and graphs to exploit the repetitive nature of textured images in order to recover fine detail 
structures. In both cases, the problem is reduced to a series of convex minimization problems that can be efficiently 
solved with a combination of variable splitting and augmented Lagrangian methods, leading to fast and easy-to-code 
algorithms. Extended experiments show a clear improvement over related state-of-the-art algorithms in the quality of 
the reconstructed images and the robustness of the proposed method to noise, different kind of images and reduced 
measurements. 

I. Formulation of the problem 

COMPRESSED sensing (CS) is founded on the principle that, through optimization, the sparsity of a signal 
can be exploited to recover it from a reduced number of measurements. This simple and yet powerful idea is 
intriguing because it seems to violate Shannon's sampling theorem. Compressed sensing is in fact the equivalent of 
Shannon's theorem from the point of view of sparsity: while Shannon states that to recover a band limited signal 
. the sampling rate must be at least twice the maximum frequency present in the signal; CS relates the sparsity of 
a signal in certain basis with the number of measurements in another basis necessary to recover it from an £\ 
minimization problem. A few definitions are necessary to understand the formulation of the CS problem. 

We say that a signal u G R n is s-sparse in the basis or dictionary ^ if it can be expressed by s non-zero 
coefficients in that basis, i.e. ||^/n||o = s; while u is compressible if most of the energy in ^u is contained in its 
largest s coefficients. Given $ and ^> two orthobasis or dictionaries of M n , the CS proplem is formulated as the 
reconstruction of a signal u G R n , sparse in basis 1 f, from m < n linear measurements / in the sensing basis <£. 
Ideally we should measure the n projections of u in basis <£, that is §u, but we only observe a small subset / = Au 
of size m < n. The sampling matrix A = results from the combination of the sensing basis <E> and the matrix 
R that extracts the corresponding to the measurements in /. Consequently, the system / = Au is undetermined 
and the sparsity of the signal u must be exploited to "invert" the problem and obtain a correct reconstruction. The 
obvious strategy would be to recover the sparsest u agreeing with the measurements, that is, to solve the following 
non-convex problem 

min ll^ullo s.t Au = f. (1) 

u 

Problem £[]) is NP-hard due to the £q norm and only approximate solutions can be used in real applications. Relaxing 
the £q norm to £\, problem (Q]) becomes the convex problem 

min s.t Au = f. (2) 

u 

One of the key results in CS CD proves that (0 exactly recovers s-sparse signals with an overwhelming probability 
when the number of measurements m > clog n (with small constant c). In addition, if the sampling matrix A verifies 
certain restricted isometry condition, then © actually recovers the signal associated to the s largest coefficients of 
u in basis i.e. exact recovery for s-sparse signals and recovery of the s-sparse £2 approximation for compressible 
signals. 
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When the measurements are contaminated with noise, the constraint Au = / on the measurements is relaxed. In 
particular, under Gaussian noise the recovery is given by 

minll^dli s.t \\Au - /|| 2 < cr, (3) 

u 

where a is related to the noise level. From optimization theory |2), we know that © is equivalent to 

min \\^u i + 7T \\Au-f \\i (4) 
u 2 

in the sense that solving any of the two determines the parameter (a,a) in the other and both have the same 
minimizer. 

Designing the sparsifying basis depends on the signal at hand. For images a common choice are orthogonal 
wavelets or the discretized total variation (TV) seminorm. TV assumes that the edges of an image are sparse and 
it is extensively used in inverse imaging problems as a regularizer. In its continuous formulation TV is a convex 
functional and its usual discretizations preserve that property. In CS H^uHi is then substituted by the regularizing 
term J(u) = \\u\\bv in an abuse of notation. 

Without loss ofgenerality in this paper we adopt the lagrangian formulation (0]), use random Fourier samples 
as measurements l_J and choose total variation as sparsity criterion; but the proposed algorithm could be equally 
applied to other basis or dictionaries as investigated in 0-|6]|. The CS recovery problem that we consider is then 

mmJ(u) + ^\\Au- f\g (5) 

With this formulation of CS recovery in imaging, we introduce an additional term in (f5]) inspired by image denoising 
techniques [7]. The resulting model exploits the geometry of the image to improve image recovery by aligning the 
normals associated to the levels sets of the image with the reconstructed signal. Our first contribution is therefore 
the introduction of a term for CS recovery based on geometric properties intrinsic to images. Our method can be 
beautifully extended to non local operators in order to recover textured images. In this case we exploit the geometry 
of the graph defined by the non local operators to recover finer details and structures of the images. This observation 
is a key contribution of our work because it can be easily adapted to improve existing non local denoising and 
deblurring methods, not only CS recovery. Finally, it is also important to mention that the proposed CS recovery 
model is based in the solution of two convex optimization problems and therefore can be efficiently solved with 
fast and easy-to-code algorithms. 

The rest of the paper is organized as follows. After formulation the problem in this section, we present our 
method in Section [TT] and explain its relation to similar techniques in Section [Till Our method is then extended to 
non local operators in Section |Tv] Section |V] presents the associated minimization problems. Finally, experiments 
are presented in Section |VI] and conclusions drawn in Section IVIII 



II. CS WITH RECOVERED NORMALS 

The main idea behind our method is that the recovered normals of an image can significantly improve the CS 
recovery results. This observation raises two main questions: how to recover normals robustly and accuratly from 
CS measurements and how to introduce the estimated normals in the CS recovery. The answer that we propose is 
a two-step iterative method. 

In the first step of each iteration, we estimate the normal vectors^ n associated to the level set curves of the image 
by solving a vectorial TV model. Once the normals are estimated, we find an image that fits the measurements, 
the estimated normals and the sparsity criterion. The process is then iterated and can be summarized as 

f n k = argmin^i^! J w {n) + ^ ||n - h\\l 

\ Uk = argmin u J(u) + 7 < div n k ,u > +^\\Au - fWl 

On the following subsections we will detail each of these two steps, which both reduce to convex optimizations that 
can be efficiently solved. Combining the two stages into one would lead to a non convex model of higher order and 

'The proposed matrix A satisfy the restricted isometry condition with high probability and is therefore a comon choice in MRI imaging 
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the resulting minimization would be slower and suffer from local minima. A two step method is computationally 
more efficient in the same way than splitting variables in Section [V] helps solving the minimization problems and 
leads to closed form solutions. The drawback of this two step procedure is the lack of rigorous theory and proof of 
convergence of the resulting algorithm. However, experimental results show that a single iteration of our method 
already improves the standard recovery (f5]), while the peak performance is attained after a few iterations. 



A. Estimation of level set normals 

Each iteration, the normals of the image are estimated in two steps. We first obtain a noisy point-wise estimate 
n from the previous solution u k ~i and we then regularize it to obtain n k . 

Given an estimate of the underlying image u k -i, the normal vectors associated to its level set curves are defined 
as 

|Vufc_i|' 

Denoising of that first estimate of the normals n is done with a combination of the vectorial ROF model [HI with 
the constraint |n| < 1. In particular we define the vector field n k = (n x ,n y ) k as the solution of the following 
variational problem. 

7 I \ I A* 1 1 - 1 1 2 , Ml - \\2 /ox 

mm J w {n x ,n y ) + -\\n x - n x \\ 2 + -\\n y - n y \\ 2 (8) 
|n|<l 2 2 y y 

where J w (n x ,n y ) is the extension of the weighted TV seminorm to vector fields and w = g(\Vu k ^i\) is an edge 
detector designed to verify w pa near the edges and w « 1 on flat regions of Uk-i- 

By weighting the TV seminorm with an edge detector w = g (\Vu k _i\), we encourage the edges of the regularized 
solution to coincide with the main edges of the noisy signal u k -i. To be robust against false edges, we use the 
robust edge detector proposed by Black, Sapiro and Marimont in where an statistical interpretation of the 
edge-stopping functions of anisotropic diffusion 1 10j is given. In this statistical interpretation, edges are considered 
outliers in the normal distribution of \Vu\ associated to noisy piece- wise constant regions and the edge-stopping 
functions g (|V«|) are derived from error norms robust to outliers. The edge detectors therefore have a parameter 
a that acts as a soft-threshold in the detection of outliers and can be estimated a priori from the values of |Vu| in 
the image. Based on the results of (9), we define 

sW-0 (l ~ f)2 !1 £ " » 

[0 \x\ > a 

with a = 1.4826 median ( | Vu — median(|Vtt|)|). 

The constraint |n| < 1 in © corresponds to a relaxation of the condition |n| = 1 inherent to the definition of 
normals. It is numerically necessary in flat regions, where Vu = and we cannot numerically normalize the 
gradient vector. 



B. Matching normals and CS measurements 

Once the normal field is computed, we find an image that matches this field by including the term — < 
rj.fc,Vu > in the standard CS recovery model (f5]). This term tries to maximize the alignment of the estimated 
normals of the signal with the normals of the reconstruction t^t. The resulting recovery model is 

u fe+ i = argmin J (u) -7 < n k+u Vu > +—\\Au- f\\\ . (10) 

u Z 

Taking into account that the divergence div is the adjoint operator of the gradient V, the previous minimization 
can be rewritten as 

u k+ i = argmin J (u) + 7 < div n k ,u > +—\\Au- f\\%. (11) 
u Z 

Our method then exploits the geometry of the image in the recovery process and obtains better regularization 
properties than standard TV. In particular the proposed model preserves edges like TV, by encouraging the gradients 
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to be sparse with J(u); but is also able to recover smooth regions by aligning the gradients of the reconstruction 
with the smoothed normals with the term < rik+i, Vu >. 

In principle we could also use a smooth estimate of the gradients v = Vitfc_i instead of n = |y^~*| to align 
the gradients of the reconstructed signal. However, discontinuities of the image would have a contribution to v 
proportional to their jump and the resulting term < v,Vu > would give different weights to discontinuities of 
different sizes. From a geometric perspective, if we want to recover the shapes of the image independently of their 
contrast we need to consider the normal vectors derived from its level sets. By the use of level sets, we treat all 
shapes equally and the term 7 < n^ +1 , Vu > only accounts for geometric quantities. 

III. Related approaches in CS 

The method that we propose is inspired by image denoising and impainting methods Q, |[TTI that align an 
estimate of the normals with the reconstructed image. In the context of image denoising, Lysaker, Osher and Tai in 
first regularize the unit gradient of the noisy image and then improve reconstruction by fitting this gradient into 
the regularized vector. The resulting method outperforms the ROF model |[T2l and similar higher order PDE methods 
lfT3ll . Dong et al. in lfl4l improve this model by regularizing the angles instead of the vectors and introducing an 
edge indicator as an extra weight. In image inpainting, an equivalent two-step method is proposed by Ballester et 
al. in [11], later improved with the divergence free constraint by Tai in lfl31 and adapted to image decomposition 
and denoising in |fT6l . In general, processing the normals to improve reconstruction has also been used in shape 
from shading [17] and mesh optimization lfl"8ll . However, this information has not been exploited before for CS 
image recovery. 

In the CS field, several methods have been proposed to improve the quality of the i\ recovery. For general 
signals, greedy algorithms |[T9l - ll2ll and t v < p < 1 minimizations ll22l . |[23l approximate the solution of the £q 
problem £[)) and improve its sparsity; but the resulting minimizations are not convex, the algorithms are slow and 
suffer from local minima. To improve the sparsity of l\ recovery © without increasing its complexity, Candes, 
Wakin and Boyd B4l proposed an iterative process solving a weigthed t\ problem at each iteration. The weights 
are defined inversely proportional to the value of the recovered signal in the previous iterate, approximating the 
behaviour of the £q norm and promoting sparser signal recovery. The resulting method efficiently solves a convex 
problem at each iteration, experimentally improves signal recovery and has been adopted for image processing with 
TV regularization in the edge-guided CS of Guo and Yin (25). Edge-guided CS incorporates information about the 
magnitude of the gradient in the recovery process and it is therefore realted to our method. However, we propose an 
additive method more robust to noise and exploit both the magnitude and directional information of the gradients. 

CS recovery of images has also been improved modifying the data term \\Au — /||| inspired by image denoising 
techniques. In particular, the Bregman iterations proposed by Osher et al. in [26] for image denoising and deblurring 
have been applied to CS in (27). He et al. in E71 use Bregman iterations to improve CS image recovery for phantom 
MRI data, but fail in the recovery of real images due to the additional difficulties of reconstructing a signal from 
partial measurements compared to the original denoising problem. For the particular case of TV regularisation, the 
fust Bregman iteration has a geometric interpretation similar to the second step of our recovery method. However, 
Bregman iterations do not include a regularization step for the normals and therefore fail for noisy and real MRI 
signals. 

In the following, we summarize each of these to methods and clarify their relationship with our technique. 
A. Edge-guided CS 

Edge-guided CS [25] improves recovery of MRI images by exploiting edge information with an iterative process 
that weights TV with an edge detector associated to the image recovered in the previous iteration. The key idea is 
that edges correspond to locations where |Vu| is large, TV corresponds to the £± norm of the norm of the gradient 
and therefore an inverse edge detector can be used to re-weight TV and approximate the £q norm in a similar fashion 
to the re-weighted l\ of Candes, Wakin and Boyd lT2"4l for general signals. The method starts with the standard 
CS solution ([5]) to obtain a first estimate of the image u\ and its edges. It then defines the weights w\ = g(\Vu\\) 
inversely proportional to |Viti|) in order to recover an image with sparser edges at the second iteration by solving 
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the re-weighted TV problem. The process is iterated, leading to the following two step algorithm: 

f u fc+ i = argmin„ J Wh (u) + § 
\ w k+1 = g{\Vu k+1 \) 

There is no stopping criterion or guarantee of convergence for this iterative process and usually, after a few iterations 
the reconstruction does not improve or even degrades. In fact, the multiplicative model of edge-guided CS is very 
sensitive to false edge detection. In particular, if an edge is detected in a wrong location, the weight associated 
to it on the next iteration will encourage an edge on this location and CS recovery will degrade with any new 
iterations. The iterative re-weighting process is designed to improve sparsity of the signal and recovery of piecewise 
constant functions, but it fails in the recovery of smooth regions in images. Compared to our method, edge-guided 
CS incorporates only information about the magnitude of Vw, while we also use its directional information; it does 
not include a regularization step for the detected edges and it is specially designed for piecewise constant images. 



B. Bregman methods 

We also share similarities with Bregman methods, whose original idea was to restore normals and image intensity 
simultaneously. However, Bregman methods cannot recover normals as accuratly and robustly as our method because 
they do not regularize the estimated normals. Our improvement is at the price of loosing global convexity. 

Bregman iterations substitute the minimization problem © for a sequence of convex optimizations substituting 
J(u) for its Bregman distance to the previous iterate. In particular, the first Bregman iteration has a geometric 
interpretation closely related to our method. Starting with u = 0, v = 0, the Bregman iterative process can be 
summarized as 

u k+ i = argmin u J(u) + §||/ + u*. - Au\\l 
v k+ i = v k + f - Au k+ i 

While their first iteration corresponds to the standard CS model (f5]), their second iteration implicitly exploits the 
normals of the image recovered at iteration one to improve the recovery. For simplicity, here we show the connection 
to our method with the continuous formulation, where A(-) is the continous functional operator of CS and A* its 
adjoint. For the first iteration u = 0, v = and the method solves 

m = argmin / |Vu| + ^\\f - A(u)\\l, (14) 

The optimality condition associated to (fT4)) derived from its the Euler-Lagrange equation is 

div = -aA*{u x ) (/ - A{ Ul )) (15) 

where n\ = j^jj correspond to the normals of u\. At the next iteration we can introduce a term < ni,Vu > 
aligning the normals of the reconstructed signal with the estimate of the normals from the previous iteration, that 
is 

U2 = argmin / iVtil— < n\, Vtt > +^11/ — A(u)||n (16) 
" Jn 2 

Integrating by parts and substituting divni in (031 ) we have 

— < rii, Vu >=< divni, u >= — < aA* (/ — A(u\)) , u >= 

-a < f - A{u\),A(u) >= -a < vi,A(u) > (17) 

with v i = / — A{i±i) as defined in (fT"3T ). If we susbstitute (fTTT ) in the minimization (fT6l ) and group together the 
terms with A{u), we end up with the Bregman update rule 

u 2 = f \v u \ + ?L\\f + Vl -A(u)\\l (18) 
Jn z 

For the rest of iterations the geometric interpretation of the update is lost. Compared to Bregman iterations, our 
method explicitly uses the normals in the recovery model for all iterations, not only the second one, and it is not 
restricted to TV regularization. Indeed, this geometric interpretation is only possible for the TV term J(u), while 
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our method can be used with any sparsifying basis. We are also more robust to noise thanks to the regularization 
step and, unlike the Bregman iteration, experimentally improve the reconstruction model © for both phantom and 
real MRI data. In addition, our method extends to non local operators to exploit graph geometry and recover fine 
details in textured images. 

IV. Extension to non local methods 

Total variation regularization is designed to recover images with sharp edges but, as other methods based on 
local gradients, it is not suited for textured images with fine structures. In this section we extend our method to 
textured images using both a non local TV regularization and a term aligning the estimated non local normals with 
the non local gradients of the reconstructed image. 

A. Non local operators 

Non local TV is a variational extension of the non local means filter proposed by Buades, Coll and Morel for image 
denoising [28 j. Non local means exploits the repetition of patterns in natural and textured images to reconstruct sharp 
edges as well as fine meaningful structures. That principle is the basis of non local regularization methods in imaging, 
which outperform the classical methods by incorporating global information in the regularization process. In |[29l 
Gilboa and Osher use graph theory to extended the classical TV to a non local functional. In the discrete setting, 
Zhou and Scholkopf |[30l and Elmoataz et al. PTI use graph Laplacians to define similar non local regularization 
operators. The resulting non local methods have been applied to image denoising 1291 , segmentation ||32l , 1331 , 
inpainting [34], deconvolution and compressive sensing |[35l . 

We adopt the discrete formulation of the continuous model presented in |2"9l . In this non local framework we 
consider the image domain as a graph G = (O, E); where Q is the set of nodes of the graph, one for each pixel in 
the image, and E is the set of edges connecting the nodes. The edge connecting nodes i and j is weighted with a 
positive symmetric weighting function w(i,j) that represents the distance between the two nodes in graph terms. 
Consequently, two pixels i and j spatially far away in the image can be considered neighbours in the graph and 
interact if w(i,j) > (we write then i ~ j). For that reason, the resulting approach is considered non local. 

Given an image u defined on the graph, the non local gradient Vqu at node i is defined as the vector of all 
directional derivatives associated to the neighbours of i, that is 

V G u = (u{j) - u(i))y/w{i,j) Vj € fl (19) 

In the graph, vectors d = d(i,j) are therefore functions defined in the domain SI x Q,. 
In this setting we define the standard L2 inner product between functions as 

< u, v >g= ^ u(i)v(i). (20) 

For vectors, we define a dot product pixel-wise 

(d-e) G (i) = ^d(i,j)e(iJ) (21) 

and an inner product on the graph 

< d,e> G = Y^{d-e) G (i) = ^ ^ d(i,j)e(i, j). (22) 

In order to have an equivalent to the TV seminorm, we define a norm function on the graph \-\g- With the previous 
definitions, the magnitude of a vector at node i is given by 

|d| G (t) = y/(d-d) G (i) = Jj2 d (iJ) 2 (23) 

The standard TV is then naturally extended to a non local version as the t\ norm of the graph norm | • \ G associated 
to the non local gradient, that is, 

TV G (u) = J G (u) = \Vgu\gW = II |VHg Hi- (24) 
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With the above inner products, the non local divergence of a vector d is defined as the adjoint of the non local 
gradient, that is 

div G d (t) = J2(d(i,j) - d(j,i))y/w(i,j). (25) 

With these definitions, if we consider only the immediate pixels as neighbours and fix their weights to w(i,j) = 1, 
then the non local TV reduces to the standard TV definition. If we consider more general neighbours by defining 
a correct weighting function like in 11281 . the non local operators incorporate global information and the standard 
regularization process is improved. The weight function therefore has an important impact in the performance of the 
non local regularizes. Inspired by 11281 . |29l , given a reference image uo we compute weighting function wo(i,j) 
measuring the difference of patches around each node as follows 

.. n-Poii)--PoU)n 2 
wo{i,j) = exp ^ , (26) 

where h is a scaling factor and Vq(i) is a patch of uo centered at pixel i. This weighting function is designed 
to reduce Gaussian noise while preserving the textures of the image. The reference image should be as close as 
possible to the true image in order to incorporate valid information related to image structures in the non local 
operators. For that reason, we initialize the weighting function in the non local methods with the standard CS 
solution ((5]) (on the following uq) and iteratively solve the non local model and update the weights with the non 
local solution. The basic non local CS recovery is then 



V Gfc 4 — estimate non local operators from u k _i 
u k = argmin u J Gk (u) + %\\ Au ~ /III 



2 (27) 



B. Proposed non local method 

Symmetrizing our local technique, we propose a two step iterative method for CS recovery. In the first step 
of each iteration, we estimate the non local normals n G associated to the level set curves of the image in the 
graph. Once the non local normals are estimated, we find an image that fits the non local normals and the CS 
measurements and iterate the process. 

In the local setting, the normal vectors associated to the level set curves of an image u are defined as n = j^jp We 
extend that definition to our non local framework and exploit the geometry of the image in the graph to improve CS 
recovery. In particular, we derive the equivalent non local normals from the non local gradient V G u by normalizing 
its components node-wise, i.e. all the components associated to node i are normalized by |V G tt| G (i). 

Given an estimate of the non local normals n G , we can include a term in the CS reconstruction (|2VT > maximizing 
the alignment of the reconstructed signal with the normals. The resulting minimization is 

u = argmin J G (u) - 7 < n G ,V G u > G +^\\Au - /||| (28) 

u Z 

Exploiting the adjoint relation of the non local divergence and gradient < n G ,V G u > G = — < div G n G ,u > G , 
we have 

u = argmin J G (u) + 7 < div G ?i,u > G +^||Au - /|||. (29) 
u 2 

As before, the process can be iterated and we obtain the following analogue to the previous two step procedure: 

V Gk < — estimate non local operators from u^-i 
div G k n Gk = argmin,, J Gk (v) + %\\v - v\ 

' ■ 2" 1 1 ^ - J II2 



u k = argmin u J Gk {u) + 7 < divG k n Gk ,u > G +%\\Au - f\\? 



with v = (l- 5 (|V G ^_i| G )) div G 

The third step of our non local method is naturally derived from our local version and the geometric interpretation 
of the non local operators. However, the regularization step of the non local normals requires careful consideration, 
as we explain next. 



s 



C. Estimation of non local normals 

The non local gradient operator, and consequently the non local normals, do not correspond to the discretization 
of standard vector fields in a grid. Indeed, V G u has a different number of components for each pixel and the 
associated direction to V G ti(i, j) depend on the relative position of the node i and its neighbour j. Therefore, 
we cannot use standard techniques to regularize these vector fields and we prefer to regularize the term div G n 
posteriorly used in the recovery algorithm. Compared to the regularization of the non local normals, we loose 
directional information, but the resulting method is simpler and faster. 

Assume that we are given an estimate of the reconstructed signal u^-i- We first compute a noisy estimate of the 
non local normals and their divergence pixel-wise and we then denoise it with standard denoising methods. In 
particular, we estimate the non local normals as 

h G = Z GUk ~\ ( 30 ) 

and compute a rough estimate of the non local divergence as 

v = (1 - 5(|V G ii fc _i| G )) div G n G , (31) 

where g(x) is a function designed to verify g w when x is large and g 1 when x is small. In fact, 5(|V G «fc_i| G ) 
acts as the equivalent edge detector presented in Section III-AI and is defined with the same expression ©. As in 
the local case, we adopt the statistical interpretation of the edge detector g (|V G u| G ) presented in |9), where the 
edges are considered as outliers in the normal distribution of |V G u| G associated to homogeneous regions. Since 
the edge detector g is derived from error norms robust to outliers, weighting our estimate of the normals with the 
function 1 — <7(|Vg«a._i|g) i n ©3 * s equivalent to soft-thresholding the non local normals when we suspect that 
the variations in Uk~i are due to noise inside homogeneous regions. 

Finally, we regularize v to obtain a smoother estimate of the non local divergence, which will be used in the second 
step of our iterative method. There are two natural approaches for this regularization: we can ignore the non local 
nature of the divergence and gradient operators and use any local model to regularize v, for instance the standard 
ROF [ 121] of equation (|32l ); or use the non local neighbours to denoise v with d33l , that is, use the natural distance 
and neighbouring relations inherent to de definition of v to denoise it. 

div Gfcn G = argmin J (v) + — \\v — v\\ 2 (32) 
v 2 

divG fc n Gfc = argmin J G (v) + -\\v - v\\ 2 (33) 
v 2 

In our experiments we obtained slightly better results with the first approach. 



V. Minimization problems 

In order to solve the minimization problems involved in each step of our method, we make use of recent advances 
in convex minimization |[36l , 11371 and apply variable splitting and augmented Lagragians ||38| to obtain efficient 
and easy-to-code algorithms. To simplify notation on this section we remove the subindexes in uu and indicating 
the iterations of our two step procedure. 

The minimizations associated to each of the local steps involve both a TV and a quadratic term similar to the 
ROF model lfl2l . Consequently, the resulting algorithms apply a similar strategy to overcome the non-linearity 
and non-differentiability of TV than the multitude of algorithms proposed for ROF. In the original ROF paper 
lfl2l . the authors derive the Euler-Lagrange PDE of the model and propose a time marching method to solve it. 
The resulting method is slow due to the constraint on the time step associated to its stability conditions. In the 
last years more efficient methods have been proposed for the ROF model due to its extensive use in imaging. A 
popular class of methods is based on the dual formulations of the ROF model, e.g. Chambolle's projection method 
||39l or primal-dual approaches ||40l - ll42"l . Other options are based on variable-splitting and equality constrained 
optimization; which is solved by quadratic-penalties ll36l . Bregman iterations 11371 . ll43l or the equivalent augmented 
Lagrangian method Il44l . In the case of CS, dual solvers are not usually adopted because they suffer from matrices 
A that are large-scale and dense. In particular for matrices corresponding to transforms with fast implementations 
(like the Fourier transform of this paper), splitting methods are a good choice because they can easily exploit fast 
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transforms to compute Au and A T u 11361 , |[37l . The algorithms that we propose fall in this last category. We rewrite 
the different problems as constraint minimizations and use augmented Lagrangians to solve them. The resulting 
Lagrangians are minimized with respect to each variable independently and the multipliers are then updated in a 
cyclic way. Since all the minimizations can be analytically solved, the resulting algorithms are extremely fast and 
easy to implement. 

Similarly, the minimization algorithms that we propose for the non local method is closely related to the 
minimization of the non local ROF model proposed in |29l , which was originally solved with a time consuming 
time marching algorithm. The non local CS problem has been solved with a combination of forward-backward 
splitting and Bregman iteration in 0311 . but for uniformity of the paper we use the same combination of splitting 
and augmented Lagrangians than in the local case to solve the non local problem 001 ). 

A. Minimizations of local normal- guided CS 

We discretize the image domain C K 2 with a regular grid of size n = n x x n y . In Q we consider images as 
scalar functions with u(i) € M and their gradients as vector-valued functions with Vu(i) G M 2 . 
We use forward differences to compute the discrete gradients and backward differences for the divergence in order 
to preserve the adjoint relationship div = — V* in the discrete setting. 
The discrete TV semi-norm is then given by 

J{u) = J2 |V«(»)I = J2 \J V *<t) 2 + VjX*) 2 (34) 

i i 

where we denote the pixel-wise norm of vectors as \d\(i) = yd x (i) + Our discretized TV is then the l\ 

norm of the function computing the pixel-wise norm of the gradient, i.e J (u) = || |Vit| 
For vector fields d = (d x ,d y ), we discretize the TV seminorm as follows 

J(d x , d y ) = J2 \/\Vd x (i)\ 2 + \Vd y (i)\ 2 . (35) 

i 

In that case we observe that it corresponds to the t\ norm of the function computing the pixel-wise norm of the 
vector of combined gradients, i.e J (d x ,d y ) = \\ \ (Vd x ,Vd y ) \ \\\. With that observation it is easy then to extend 
it to a weighted TV norm as J w (d x ,d y ) = || \W {Vd x ,Vd y ) \ ||i, where W is the diagonal matrix of weights. 

In the vector notation used in CS, we can efficiently compute the spatial derivatives multiplying the discrete 
functions arranged as a column vector with the sparse finite difference matrices V x u = D x u, V y u = D y u. 
Similarly, the discretization of the L2 inner product in SI corresponds to the usual dot product of vectors, i.e. 
< V , u >= v T u. 

1) Estimate u from CS measurements and normals: To reconstruct the image we need to solve the following 
convex minimization problem: 

min II \Vu\ ||i + jv T u + — \\Au — /|| 2 with v = divn. (36) 
u 2 

We propose an iterative algorithm to solve (l36l ) based on splitting and constraint minimization techniques. The main 
idea is to split the original problem into sub-optimization problems which are well-known and easy to solve, and 
combine them together using an augmented Lagrangian. The proposed algorithm is guaranteed to converge thanks 
to the convexity of (l36l) . 

Let us consider the following constrained minimization problem, which is equivalent to (ITTb : 

min = || |d| ||i + v T u + —\\Au - f\\\ s.t. d = Vu (37) 

u,d 2 

The proposed splitting approach makes the original problem (ITTb easier to solve because (l37l) decouples the t\ 
norm and the gradient operator V. 

Next, we reformulate this constrained minimization problem as an unconstrained optimization task. This can be 
done with an augmented Lagrangian approach, which translates the constraints into pairs of Lagrangian multiplier 
and penalty terms. Let us define the augmented Lagrangian energy associated to (l37l ): 
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A(«,d,A) = || \d\ \\ 1 + v T u + -\\Au- f\\ 2 2 + \l(d x - D x u) 

+ \y(dy - DyU) + ^Wd.^ - D X U\\l + r ~\\dy - DyU \ \ \ (38) 

where A = (X x , X y ) are the Lagrange multipliers and r is a positive constant. 

The constraint minimization problem (137T ) reduces to finding the saddle-point of the augmented Lagrangian energy 
C\. The solution to the saddle point problem (|38l ) can be approximated iteratively by the following algorithm: 
initialize the variables and Lagrange multipliers to zero; at each iteration find an approximate minimizer of 
L\ (u, d, A fc ~ x ) with respect to the variables u, d and update the Lagrange multipliers with the residuals associated 
to the constraints; stop the process when u remains fix. As the Lagrangian L\ is convex with respect to u, d, we 
can find a minimizer by iteratively alternating the minimization with respect to each variable. The resulting method 
is equivalent to the alternating direction method of multipliers. The iterative method is summarized in Algorithm Q] 



Algorithm 1 Augmented Lagrangian method to solve (1371 ). estimating u from CS measurements and normal 
matching 

1: Initialize u, d, A 

2: For each iteration 1 = 1,2..., find an approximate minimizer of C\ with respect to variables (u, d) with fixed 
Lagrange multipliers A': 

u l = argmin£i (u, d 1 , A i_1 ^ solved in in Fourier domain (39) 
d l = argmm£i , d, A' -1 ^ solved by shrinkage (40) 

3: Update Lagrange multipliers 

=X l ~ 1 +r(4 -D x u l ) 
X L y =\ l - 1 +r(d y -D y u l ) 

4: Stop the iterative process when ^" jj^nj — - < (■ 

The next step is to determine the solutions of the two sub-minimization problems (|39l),(l40b. which can be 
computed efficiently. 

The sub-minimization problem (|39l ) can be written as follows: 

min v T u + %\\Au - /||1 + ^-\\d x + -A* - D x u\\\+ 
u 2 2 r 

r -\\d y +^-Xy-DyU-\\l (41) 

We see that it reduces to a quadratic minimization, with positive semi- definite Hessian H = aF T R T RF+r(D^ D x + 
DyDy). The optimality conditions read 

Hu = b with b = aF T Rf + D T X (rd x + Ax) + {rd y + X y ) . (42) 

Actually as R is a row selector, R T R is a sparse diagonal matrix with non-zero entries on the selected Fourier 
coefficients and we cannot assure the invertibility of H. We find an approximate solution defining the positive 
definite matrix H € = H + el n with small e > and solving the approximate system 

H t u = b + eu, (43) 

where we use the value of u from the previous augmented Lagrangian iteration to estimate u = u 1 ^ 1 . In the 
resulting system, H t is block circulant and we can use the Fourier transform to decompose it as H e = T^CF, 
with C = aR T R + rT {p^D x + DyD y ) T T + el n a diagonal matrix. Consequently, the system (l43l 
can easily be solved in the Fourier domain inverting the diagonal matrix C. In practice we use the FFT transform 
instead of doing the matrix multiplications with T and T T , which gives us a solution of complexity 0{n\ogn). 
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The minimization problem w.r.t. d corresponds to an l\ - £2 norm and can be solved by shrinkage. We first note 
that (1401 is equivalent to 

r 1 
min + ~Y^ l d W ~ z ^)! 2 ' with z = -X- Vit. (44) 

I I 

The minimization of (l44l can be done pixel-wise and the solution is given by the shrinkage operator S (z,l/r). 

d(i) = max \\z(i)\--,o\^l T i = l,...,n (45) 
I r J \z(i)\ 

2 ) Regularization of normals: To regularize the normals at each iteration we have to solve 

min || \W (Vn x , Vn y ) | ||i + —\\n x - n x \\l + ^||% - h y \\l, (46) 
|n|<l Z Z 

where W is a diagonal matrix with weights associated to weighted TV seminorm. We use the same combination 
of splitting and augmented Lagrangian techniques than in Section IV-A1I To avoid repetition, on the following we 
will simply write the form of the constraint minimization problem, the augmented Lagrangian and each of the 
subminimizations for a self-contained paper. 
Equivalent constraint minimization is 

min || \W(d,e) | ||i + -\\n x -n^ll + -11% -%|||, (47) 

n=m,\m\<l Z Z 

d=Vn x ,e=Vn„ 



with associated augmented Lagrangian 

i.e. X.v.t) = II \W(d.e) I Hi 4- 



£2 (n,m,d,e, A, = || |W(eZ,e)| ||i + — ||n — ri.|||+ 



X x (d x - D x n x ) + \y(d y - D y n x ) + -^\\d x - D x n x \\l+ 

y \\dy - DyUyWl + ~ D X Uy) + Vy^y ~ DyUy) 

~l~ 1 1 e ^ _ D X Tly 1 1 2 H — || &y — DyUy 1 1 2 + £, x (jl X ~ TH X ) 

+-^\\ n x ~ m x \\ 2 + £ y (n y - m y ) + —\\n y - m y \\ 2 . (48) 

The resulting minimization method is presented in Algorithm |2] 

The sub-minimization problem with respect to n x can be written as follows: 

^ II - II 2 1 t-T I \ 1 r m 11 11 2 

mm -\\n x - n x \\ 2 + £ x (n x - m x ) + — \\n x - m x \\ 2 

n m Z Z 
V ^ d 

+ dl<4 H K - D x n x \\l + —\\d y H \ y - D y n x - \\l (49) 

& T"d " r d 

We see that it reduces to a quadratic minimization, with positive definite Hessian H = (fJ, + r m )I n + r d D^D x + 
r^D^Dy. The optimality conditions read 

Hn x = [ih x + r m m x + L>J (r d d x + \ x ) + Dy (r d d y + X y ) - £ x . (50) 

As before, the resulting H is block circulant and we can use the Fourier transform to decompose it as H = T T CT, 
with C = (fi + r m )I n + r d T (D X D X + DyD y ) T T a diagonal matrix. We solve the linear system in the 
Fourier domain efficiently with the FFT transform. Observe that the minimization problem with respect to n y has 
the same form and can be solved with the same technique. 

The minimization problem w.r.t. d corresponds to the £\ - £2 problem 

min^ \w(i)d(i)\ + \d(i) - z(i)\ 2 , with z = —A - Vn x , (51) 

% % 

which is equivalent to the following minimization (note w(i) > 0) 

min \d(i)\ + t^tt £ \d(i) - z{i)\\ with z = —A — Vn x . (52) 
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Algorithm 2 Augmented Lagrangian method to solve (l47l) . regularizing normals 



1: Initialize n, m, d, e, A, £ 

2: For each iteration Z = 1, 2 . . ., find an approximate minimizer of £2 with respect to variables (n x , n y ,m, d, e) 
with fixed Lagrange multipliers : 



3: Update Lagrange multipliers 



n l = aigmiiiC2(n,rn l 1 ,d 

TO 


-\e l - 




-1 


v l 


-\e 


- 1 ) 


m 1 = argmin£(n i , m, <i 

m 


-\e l - 




-1 


v l 


-\e 


- 1 ) 


d l = argmin£(n', m l 

d 


,d, e l - 




-1 


v l 


-\e 


- 1 ) 


e l = argmin£(n i , m\ d l , 


e,A< 


-1 


v l 


-\t 


- 1 ) 



X 1 ^^ 1 +r d {d l x -D x n l x ) 
X l y =X l i 1 + r r (d l y - D y n l x ) 
u l x =v l ~ x + r e (4 - D x n!y) 
v l y =v 1 - 1 +r e {e l y -D y n l y ) 

A _A-\ 



C =e- 1 +r m (n-m) 



4: Stop the iterative process when ^ n n£k — - < e. 



A similar problem has already been solved in Section IV- A 1 1 with the shrinkage operator, which is now adapted to 
include the weights w. The solution is then 

d(i) = max \\z(i)\ - ^-,o\-^- i = l,...,n (53) 
L r ) \z[i)\ 

Due to the symmetry of the problems, the same minimization technique is used for e. 
Finally, the minimization problem w.r.t. m reads 

min — \^ \m(i) — z(i)\ 2 , with z = n + — ^ (54) 

|m(i)|=i 2 ^ r m 

and can be solved pixel-wise. For each pixel we have the following 2-D problem: given a point in space with 
coordinates z(i) € M? we want to find the point in the unit ball minimizing its distance to z(i). It is clear that the 
solution corresponds to the projection of the unconstrained minimizer z(i) into the unit ball, i.e 



z(i) \z(i)\ < 1 
otherwise 



m (i) = { z(i) _ it .__., ; „_ • (55) 



B. Minimizations of non local normal-guided CS 

In the discrete setting, the NL gradient is a linear operator. Arranging the image as a column vectors, it can be 
computed as a sparse matrix multiplication Vqu = Du. The matrix D € M Arxn (N = \E\ indicates the number 
of nodes in the graph) is derived from the weights associated to the edges and is usually sparse. Consequently 
d = Du G W N is also a vector column, with as many components associated to node i as neighbours this node 
has. With the vector notation, the inner product of two vectors fields d, e defined in G is then computed as 
< d,e >g= d T e. As in the continuous setting, the NL divergence div^ is derived from its adjoint relation with 
the NL gradient Vq = — div^ and, consequently, in matrix notation it corresponds to div^ d = —D T d. 

Since the minimization associated to (l32l has already been explained for the vectorial case, in the next paragraphs 
we focus on the minimizations associated to non local operators (l30l) and (l33l) . 
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1) Minimization associated to CS reconstruction matching non local normals: With the previous notation, the 
minimization problem (l30l reads 

u = arg min || |Dm|g ||i + jv T u + — \\Au — ff (56) 
u 2 

with v = dive no- This minimization is also reformulated as a constraint minimization and solved efficiently 
with augmented Lagrangians. Compared to the local minimizations, in the splitting step we require an additional 
variable s to have efficient and analytic solutions for the posterior subminimization problems. The resulting constraint 
minimization formulation of (l56l) is 



mm = || \d\ G \\ 1 +v T u + -\\Au- ff 2 s.t. <^ _ (57) 

u,s,a Z IS — U 



The Lagrangian in that case reads 

£ 3 {u, s, d, X d , X u ) = || \d\ G ||i + v T u + -AAu - ff 

+\ d T (d-Du) + T A\d- Du\\ 2 + \l{u-s) + T -^\\u - sf. (58) 

The resulting minimization method is presented in Algorithm [3) where we have also hinted the solution to each 
of the subminimization problems. 



Algorithm 3 Augmented Lagrangian method for CS reconstruction matching normals by 



1: Initialize u, s, d, Xd, X u 

2: For each iteration I = 1,2..., find an approximate minimizer of £3 with respect to variables (n, s,d) with 
fixed Lagrange multipliers A^', A^: 

u = argmin/^ti, s 1 ^ 1 , d'" 1 , X^ 1 , A^ -1 ) solved with conjugate gradients 

u 

s = arg min £3 (u l , s, d 1 " 1 , X^ 11 , A^ 1 ) solved in Fourier domain 

s 

d = arg min £ 3 (u L , s l ,d, Ad' -1 , A^ 1 ) solved by non local shrinkage 

d 

3: Update Lagrange multipliers 

X d l =X d 1 - 1 +r d (d l -Du 1 ) 
Xl=X l u - 1 +r u (u l -s l ) 

II ( —y 1 - 1 1 

4: Stop the iterative process when - — — - < e. 



The minimization w.r.t u corresponds to the following quadratic positive definite problem 

min v u + Xd (d — Du) H \\d — Du\\ + X u (u — s) H \\u — s\\ . 

u 2 2 

We find its minimizer solving its optimality conditions, which provide the following system of linear equations 

(r u I + r d D T D) u = -jv - X u + r u s + D T (X d + r d d) . (59) 

Matrix K = r u + r d D T D is sparse, symmetric and positive definite and we have efficient algorithms to invert it. 
We choose an iterative method to invert the matrix, initializing it from the previous solution to the minimization 
problem u l ~ l . In particular we use the conjugate gradient method to exploit the symmetry and positive definition 
of K, with preconditioning matrix given by its incomplete Cholesky factorization. The resulting method is very 
fast, converging to enough precision with 2 — 3 iterations of the conjugate gradient method. 

The minimization w.r.t s is also a quadratic problem which can be efficiently solved, in that case in the Fourier 
domain. The problem reads 

min %\\Au - /|| 2 + Xl(u - s) + ^||u - sf. (60) 
s 2 2 
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The optimality conditions in that case are 

(aA T A + r u I n ) s = aA T f + X u + r u u. (61) 
As before, the matrix A T A + r u I n = F T CF is block-circulant and the resulting system is diagonal in the Fourier 



domain with C = R T R + r u I n . Therefore (I6TT) can be efficiently solved with the FFT 
The minimization with respect to d is equivalent to 

min = II \d\c ||i H — — z|| 2 with z = Du — — (62) 
d 2 r d 

As in the local case, this minimization is decoupled for each pixel i as follows 



= j£*(M) + £?(<i,i)-*(<,i)) a (63) 



mm 



d*(i,j) = max{|z| G (i) - l.ojf^-. (64) 



and can be solved by a straight-forward extension of the shrinkage operator to the graph. That is, for each node 
neighbour to i the solution is given by 

2J Minimization associated to the regularization of normals: With the previous notation, the minimization 
problem (l33l reads 

v = argmin LDv g iH — \\ v ~ v 2 (65) 
« 2 



As in the local case, we decouple the l\ and £2 problems defining an additional variable d = Du and rewrite (1651 ) 
as the following constraint minimization problem 

fi 2 

min = II \d\Q ||i H — ||w — €> 1 1 2 s -t- d = Du (66) 

with associated augmented Lagrangian 

£ 4 (n, d, A d )= || |d|e ||i + - ||u - v\\l 

+X d T {d- Dv) + 7 j\\d- Dv\\l (67) 

To minimize the Lagrangian £3 with respect to u, d, we alternate the direction of minimization with respect 
to each variable and proceed as indicated by Algorithm |3l where we have also hinted the solution to each of the 
subminimization problems. 

Algorithm 4 Augmented Lagrangian method to regularize non local divergence of normals from d66l) 
1: Initialize u, d, Xd 

2: For each iteration Z = 1,2..., find an approximate minimizer of £4 with respect to variables (u, d) with fixed 
Lagrange multipliers A,/: 

v = argmin£4(w, d' 1 , A^' 1 ) solved with conjugate gradient 

V 

d = argmin£4(u i , d, Ad' -1 ) solved by non local shrinkage 

d 

3: Update Lagrange multipliers 

A d ( =A d *- 1 + r d (d l - Du 1 ) 
4: Stop the iterative process when ^ v jj~j~y — - < e. 

The minimization w.r.t v corresponds to the following quadratic positive definite problem 

min -\\v - v\\l + \ d T (d - Dv) + — lid - Dv\\ 2 2 . (68) 
v 2 2 
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We find its minimizer by solving its optimality conditions, which provide the following system of linear equations 

(fil + r d D T D) v = fiv + D T (A d + r d d) . (69) 

We find the same form of matrix K = + r^D 1 D than in d59l and, therefore, we solve with linear system d59l ) 
with the same conjugate gradient method. 

The minimization with respect to d is equivalent to d62l ) changing u for v, in particular we have 

min = d G hH d — z with z = Du (70) 

d 2 r d 

and is solved with the same adaptation of the shrinkage operator to the graph. For each node neighbour to i, the 
solution is given by 

d*(i,j) = max\\z\ G (i) - -,0}P^- V (71) 
L r d J \z\G\i) 

VI. Numerical results 

In this section we present some of the numerical results obtained with our method and compare it to other 
techniques. We compare the local version of our method to standard CS recovery algorithm (HJ) and to the edge- 
guided CS proposed in |[25l . The non local version of our model is compared to the non local CS recovery (l27l ). 
which does not take into account the geometric information of the non local gradients into the recovery process. 
For the non local case, in our model we have regularized the divergence of the normal with the standard ROF 
model. 

We use partial Fourier measurements for our reconstruction and perform radial sampling on R with different 
number of measurements in relation to the size of the signal (we specify it with the ratio — ). For a fair comparison, 
we have used the same robust edge detector (O for the edge-guided CS and our method and we have implemented the 
minimizations with the same splitting and augmented Lagrangian techniques for all the methods. The parameter a, 
which is related to the noise present in the CS measurements, has been manually tuned to obtain best reconstruction 
with the standard CS recovery models © and (|27T ) and used with the other methods. The other parameters of our 
model 7, fi have also been chosen manually to obtain good CS recovery in terms of SNR. We have observed that 7 
(which controls the weight given to the alignment of the normals) takes similar values for the same kind of images 
(textured or brain IRM images) and remains stable for different sparsity and noise levels. On the other hand, the 
parameter fi controlling the smoothness of the estimated normals decreases when the number of measurements 
decreases or the noise level increases because the partial reconstructions and the estimated normals are noisier and 
require more regularization. 

In a first set of experiments we test our method with MRI images, first with the Shepp-Logan phantom from 
Figure Q] and then with a real MRI brain image from Figure |2l where we include also the results with back- 
projection, i.e. the solution to / = Au with smallest I2 norm. 

Table U show the quantitative results of the different CS reconstruction methods for MRI images. Our method 
always outperforms the standard TV reconstruction and the edge-guided CS technique. In the experiments, both 
the edge-guided CS and our proposed method are initialized with the TV solution and, therefore, always improve 
its reconstruction. Comparing the gains of these two methods with respect to the TV reconstruction, we observe 
that our method more than doubles the gain of edge-guided CS. Figurefl]|3] show qualitatively the improvement 
of our method over TV reconstruction. In the case of the phantom we are able to better reconstruct the phantom 
with fewer measurements both in the local and non local case, while with a real MRI image our reconstruction is 
able to capture better non-dominant edges of the white-grey matter interface. In Figure [3] we explicitly compare 
the normals associated to the TV solution and the regularized normals of our local reconstruction for the real brain 
MRI image. We observe that our method is able to better reconstruct the normals, and therefore the shapes, of the 
image (to clearly see the differences please check the PDF, not the print out). 

Performance improves with non local regularization, with our method outperforming the non local CS reconstruction 
for all the experiments. As expected, the gain of our method compared to TV is lower than in the local approach 
because we loose part of the directional information of the normals by denoising their divergence instead of the 
vector fields. For each image we also added different levels of Gaussian noise (a n ) to the signal to investigate the 
robustness of our method to noise. Results are shown in table JI] We observe that we are more robust to noise than 
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Image 


m 
n 


local CS 
TV edge CS normal CS 


non-local CS 
TV normal CS 


Phantom 
Phantom 
Brain 
Brain 


8% 
12% 
12% 
20% 


7.33 dB 7.37 dB 12.78 dB 
38.60 dB 45.33 dB 56.14 dB 
17.14 dB 17.38 dB 17.71 dB 
22.16 dB 22.35 dB 23.82 dB 


28.28 dB 33.13 dB 
61.84 dB 74.57 dB 
18.96 dB 20.39 dB 
23.13 dB 24.12 dB 



TABLE I: Comparison of CS reconstruction for MRI images. The first three columns show the results with the 
standard TV term in the regularization: TV stands for the model (f5]), edge CS for (fl2l ) and normal CS for our 
method Q. The last two columns correspond to the definition of NL TV: NL-TV corresponds to the standard non 
local CS recovery (|27T ) and NL normal CS for the proposed non local method (l30l ). 

edge-guided CS (which in fact does not improve the TV reconstruction for noise levels a n = 15%, a n = 10%) 
thanks to regularization step on the estimation of the normals. As before, non local regularization improves CS 
reconstruction, we observe that our non local method outperforms again the non local TV and is also robust to 
noise. 



Image 

S = 12% 

n 


noise 

On 


local CS 
TV edge CS normal CS 


non-local CS 
TV normal CS 


Phantom 


5% 
10% 
15% 


11.90 dB 11.91 dB 12.90 dB 
8.37 dB 8.38 dB 9.44 dB 
6.59 dB 6.59 dB 7.28 dB 


17.92 dB 18.36 dB 
12.15 dB 13.03 dB 
10.09 dB 10.27 dB 


Brain 


5% 
10% 
15% 


13.37 dB 13.36 dB 13.78 dB 
10.88 dB 10.88 dB 11.57 dB 
9.89 dB 9.89 dB 10.48 dB 


14.86 dB 15.00 dB 
12.31 dB 12.50 dB 
10.94 dB 11.19 dB 



TABLE II: Comparison of CS reconstruction for noisy MRI images with 12% of samples and different levels a n 
of gaussian noise. The first three columns show the results with the standard TV term in the regularization: TV 
stands for the model (HJ, edge CS for (f72l and normal CS for our method ©. The last two columns correspond to 
the definition of NL TV: NL-TV corresponds to the standard non local CS recovery (|27T ) and NL normal CS for 
the proposed non local method d30"l) . 

In a second set of experiments, we tested our method with natural images containing textures, where edge 
detection by itself is a difficult task and the images can not be considered piecewise constant. With these images, 
the local regularization looses all texture information, while the non-local approaches can recover repetitive patterns 
and better exploit the geometrical information of the image. Results with our method are presented in table [[III 
with some of the reconstructed images shown in figures |4]|9] 

A quantitative comparison of the different methods with textured images is presented in table |III] We observe 
that the inclusion of an edge detector in edge-guided CS does not improve the TV reconstruction because the 
partially reconstructed images are not accurate enough to detect edges and the weighted TV term of edge-guided 
CS encourages edges in wrong positions. That effect is not observed in our method because it is additive and not 
multiplicative and it exploits the directional information of the regularized normals, which can partially capture 
texture information better than an edge detector. As a consequence, our local method always outperforms the TV 
reconstruction and edge-guided CS methods. For the non local regularizations our method outperforms non local 
TV, but the gain in some cases is negligible (fingerprint and baboon images for a ratio of measurements — = 12% 
or 20%). In fact, the non local methods require a good estimate of the reconstruction to initialize the non local 
gradient and divergence operators. Since our method requires both gradient and divergence to estimate the non local 
normals and align them with the reconstruction, we can only improve the non local TV reconstruction when the 
initialization (in our case we use the standard TV solution) has a minimum level of accuracy. The fact that more 
measurements are required for the fingerprint of baboon images is coherent with CS theory, as these images have 
finer details and are less sparse than Lena or Barbara in terms of total variation. In figure [4][5] we can qualitatively 
observe the advantages of our method in comparison to local and non-local TV reconstruction for textured images. 
In the local case we avoid the staircase effect, which is clearly visible in the TV reconstruction of Lena's cheeck. 
In the non local case, we also capture better slowly varying textures changes, see for instance the different shadows 
in Lena's skin or hat. In both cases this improvement is due to the regularization of the level set normals of the 
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image, which we exploit for CS recovery with our two step procedure. 



m 
??. 


image 


local CS 
TV edge CS normal CS 


non-local CS 
TV normal CS 


12% 


Lena 
Barbara 
fingerprint 
baboon 


14.53 dB 14.47 dB 14.86 dB 
13.35 dB 13.31 dB 13.59 dB 
4.13 dB 4.11 dB 4.13 dB 
7.40 dB 7.25 dB 7.40 dB 


15.82 dB 16.79 dB 
15.00 dB 15.52 dB 
5.97 dB 5.98 dB 
7.65 dB 7.65 dB 


20% 


Lena 
Barbara 
fingerprint 
baboon 


18.44 dB 18.36 dB 19.27 dB 
16.71 dB 16.62 dB 17.13 dB 
5.70 dB 5.62 dB 5.70 dB 
9.13 dB 8.91 dB 9.14 dB 


19.95 dB 21.09 dB 
18.37 dB 18.93 dB 
9.03 dB 9.07 dB 
9.63 dB 9.74 dB 


39% 


Lena 
Barbara 
fingerprint 
baboon 


25.39 dB 25.30 dB 26.71 dB 
20.83 dB 20.68 dB 21.36 dB 
12.02 dB 11.84 dB 12.03 dB 
13.30 dB 13.14 dB 13.41 dB 


26.39 dB 27.51 dB 
24.68 dB 25.33 dB 
14.52 dB 14.56 dB 
13.44 dB 13.82 dB 



TABLE III: Comparison of CS reconstruction for textured images. The first three columns show the results with 
the standard TV term in the regularization: TV stands for the model (f5]), edge CS for (fT2l) and normal CS for our 
method The last two columns correspond to the definition of NL TV: NL-TV corresponds to the standard non 
local CS recovery (TZTT ) and NL normal CS for the proposed non local method (l30l ). 



VII. Conclusions 

We propose a normal guided compressed sensing recovery method to recover images of higher qualities from fewer 
measurements. The normal vectors of image level curves have been exploited in denoising and inpainting algorithms, 
but in compressed sensing this information is embedded in the measurements and state-of-the-art recovery algorithms 
have just neglected it. To extract this geometric information we alternatively estimate the normals of the image level 
set curves and then improve the compressed sensing reconstruction matching the estimated normals, the compressed 
sensing measurements and the sparsity constraints. With our approach, the geometric information of level contours 
are incorporated in the image recovery process. The proposed method is also extended to non local operators to 
recover textured images and could also be applied to improve existing non local denoising and deblurring methods. 
Our numerical experiments show that the proposed method improves image recovery in several ways: it is able to 
recover sharp edges as well as smoothly varying image regions, avoiding the staircase effect in the case of total 
variation reegularization; it is robust to noise and the sparsity of the signal and relies on efficient minimization 
techniques to obtain a fast and easy-to-code algorithm. 
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(c) TV reconstruction, 7.33 dB (d) non local TV reconstruction, 28.28 dB 




(e) proposed local method, 12.78 dB (f) proposed non local method, 31.26 dB 

Fig. 1: Reconstruction of Shepp-Logan phantom from 8% of measurements in Fourier domain. 



(a) Original image 



(b) Back-projection, 9.25 dB 




(e) proposed local method, 18.56 dB (f) proposed non local method, 20.39 dB 

Fig. 2: Reconstruction of MRI brain image from 12% of measurements in Fourier domain. 
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(a) Original image (b) Proposed local technique 



Fig. 3: Zoom on reconstructed brain MRI image from 12% of measurements in Fourier domain. We superpose the 
reconstructed signals with the normals associated to their level sets for the standard TV solution (left) and for the 
local version of our method (right). Our method is able to better reconstruct the normals and shapes of the image. 



(a) Original image 



(b) Back-projection, 8.30 dB 




(c) proposed local method, 13.59 dB (d) proposed non local method, 15.52 dB 

Fig. 6: Reconstruction of Barbara from 12% of measurements in Fourier domain. 
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(c) proposed local method, 9.14 dB (d) proposed non local method, 9.74 dB 

Fig. 8: Reconstruction of baboon from 20% of measurements in Fourier domain. 



(a) Original image 



(b) Back-projection, 4.69 dB 




(c) proposed local method, 5.70 dB (d) proposed non local method, 9.07 dB 



Fig. 9: Reconstruction of fingerprint from 20% of measurements in Fourier domain. 



